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Starting from the Random Phase Approximation (RPA), we generalize the schematic model of 
separable interaction defining subspaces of ph excitations with different coupling constants between 
them. This ansatz simplifies the RPA eigenvalue problem to a finite, small dimensional system of 
equations which reduces the numerical effort. Associated dispersion relation and the normalization 
condition are derived for the new defined unknowns of the system. In contrast with the standard 
separable approach, the present formalism is able to describe more than one collective excitation 
even in the degenerate limit. The theoretical framework is applied to neutral and singly charged 
spherical sodium clusters and Ceo fullerene with results in good agreement with full RPA calculations 
and experimental data. 

PACS numbers: 


I. INTRODUCTION 

Recently 0,0 we have investigated the collectivity of 
the Giant Dipole Resonance and of the low-lying modes 
(Pigmy Dipole Resonance) in nuclear systems using the 
separable RPA. In contrast with the standard schematic 
model 0 which employs separable residual interaction 
with a single specific coupling constant between nucle¬ 
ons, we have assumed the existence of two subsystems 
of particle-hole states and, consequently, three distinct 
coupling constants. 

In the present work we generalize this approach to a 
variable number of subsystems of particle-hole pairs, al¬ 
lowing us to describe an even higher number of collective 
resonances - a feature which is relevant to systems such 
as metallic clusters. 

The optical response of clusters has been intensively 
studied in the past three decades from both an experi¬ 
mental and a theoretical point of view (for reviews see 

Si). 

The most pregnant effect in the optical response of 
large (metal) clusters, is the surface plasmon. This rep¬ 
resents a dipolar motion of the electron cloud against the 
ionic background, similar with the Goldhaber-Teller [0 
description of the Giant Dipole Resonance in nuclei. It 
is interpreted as a collective phenomena given by the su¬ 
perposition of single particle effects at nearly the same 
frequency. Another collective resonance is the so called 
volume plasmon which, from a macroscopic view is char¬ 
acterized by density variations inside the cluster volume 
as in the Steinwedel-Jensen model 0. 

These are gross properties of the optical spectra, ex¬ 
hausting large fractions from the total sum rule in large or 
spherical shaped clusters. Still, spanning different sizes, 
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the shell effect becomes stringent [8j, similar with the 
atomic case, allowing for the existence of many deformed 
clusters in between the spherical ones. Moreover, for clus¬ 
ters with a small number of atoms, the details of the ionic 
background become important. For all that reasons, we 
face in optical spectra relevant single particle excitations 
or fragmentations of the collective surface plasmon. The 
latter is connected with the phenomenon of Landau frag¬ 
mentation. 

Unlike in nuclear physics, where Hartree-Fock (HF) is 
the main microscopic theory, the dynamics of electrons in 
clusters is usually investigated by means of the Density 
Functional Theory (DFT) 0-[l(| which is implemented 
through the Time-Dependent Kohn-Sham equations Q 
The extension to the time-dependent case describes the 
properties of excited states and its small amplitude limit 
represents the random phase approximation (RPA) cor¬ 
respondent of the HF derived RPA. 

Even though RPA can be easier to implement than 
full HF or DFT, it still poses a difficult numerical task. 
This is reflected in the dimension of the space of single 
particle excitations, given the fact that for each pair of 
ph states, two 6 D integrals should be computed. The 
difficulty arises especially in deformed clusters, where the 
space of the single particle states lacks degeneracy and 
the spherical symmetric is removed. 

To overcome such challenges, there has been an ac¬ 
tive research in the past decades in the art of simplified 
RPA methods. The lowest level of approximation is given 
in the frame of sum rule techniques fl2j where a single 
resonance is described. A more pertinent description is 
achieved with the local RPA method [10 where a set of 
local operators is used to describe an equal set of col¬ 
lective resonances, but the Landau damping features are 
lost. The original separable ansatz [3] is used in cluster 
physics as the vibrating potential model EMI- This has 
the main disadvantage that only one collective excitation 
is obtained and the natural fragmentation of the surface 
plasmon is not properly described. 
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A more refined method is provided by the self- 
consistent separable RPA (SRPA) [H ill which com¬ 
bines the idea of more local operators and the separabil¬ 
ity ansatz of residual interaction. 

In this work we generalize the separability of the resid¬ 
ual interaction in the frame of DFT-LDA using a single 
operator but many coupling constants prescribed by N 
blocks of ph excitations, method which will be referred to 
as N Block Separable RPA (./V-BSRPA). Our formalism 
has many similar points with the one developed in 0 , 
but, as we shall see, the interpretations are fundamen¬ 
tally distinct. While SRPA describes a system of coupled 
collective motions, Af-BSRPA accounts for a collection of 
coupled subsystems of oscillators. 

The paper is organized as follows: In Section [TT] a short 
derivation of DFT derived RPA is given and the formal¬ 
ism of ./V-BSRPA is developed. In Section Hill the analyt¬ 
ical results are applied for some sodium clusters and Ceo 
fullerene. 


II. THEORY 
A. RPA from DFT 

The Time-Dependent DFT 0 within the adiabatic lo¬ 
cal density approximation (LDA) 0 can be represented 
through the Kohn-Sham equations, which for a system of 
electrons described by single particle orbitals {<^(r, t)}, 
read: 


h[p(r,t)\<t>j( r,t) =iH3 t ^j(r,t), (1) 

with 

p( r 1 *)=Efe'( r ’ t )i 2 ' ( 2 ) 

jEocc 

j e occ indexes the occupied single particle states, p( r, t) 
is the total density of the particles and h[p{ r, t)] is a mean 
field Hamiltonian incorporating, besides the kinetic term, 
the Coulomb interaction between electrons, a local poten¬ 
tial for the exchange-correlation effects and an external 
potential. Since we are interested in the normal modes 
of the electron system, we consider that the external po¬ 
tential is being given only by the Coulomb interaction 
between electrons and ions. 

In the stationary case, the KS equations become eigen¬ 
value problems h[po(r)]ipj(r) = £j<pj(r). Each orbital 
can be expanded around its stationary value 4>j (r, t) = 
(ifj (r )+5(f>j (r, t))e~ lSjt / h and with this expansion we can 
linearize the Hamiltonian and the density as: 


Mm) = 2 E # e kj( r )M(M)], 

jdzOCC 


h[p] = h 0 [p 0 \ + 


Sh 

Sp 


Sp. 

Po 


(3) 


Since we are interested in the oscillating behaviour of 
5<f>j we search for the harmonic coefficients of the ex¬ 
pansion in the basis of both occupied and unoccupied 
single-particle states: 


M(M) = E M*) (x jlie -^ + Y*^) . (4) 

Replacing the linearised forms [3] and the expansion |4] 
of 5(f)j{v,t) in eq. [T] and integrating the equation after 
multiplication with </?„, we obtain the usual form of the 
RPA eigenvalue equations 


£ iX[ k) + E( A b X i fe) + BijY^ 

3 

stY^ + Y / (B* j X^ k) +A* j Y} k) 

3 

with J2j(\Xj k) \ 2 — I^El 2 ) = 1 as normalization con¬ 
dition . In deriving the above equations, we used the 
compact notations Xj^ = Xf, Yj^ = Yi, — £j = £j, 
Aij = (tpvPj \Sh/Sp\pi<p^) and Bij = (ip^ip^Sh/Splpiipj). 
The superscript (k) labels the solutions of the eigen¬ 
value problem. For legibility we shall suppress this nota¬ 
tion whenever possible. In literature, the results of this 
derivation are also known as the linearized TDDFT. More 
detailed discussions and related derivations can be found 



It might seem that, unlike in the HF derived RPA, the 
system of eqs [5] describes transitions from occupied states 
into other occupied states, fact which violates the Pauli 
principle, but in practice, since these terms cancel each 
other, we will set them to 0. While in the standard RPA 
the residual two-body interaction is written as a differ¬ 
ence between direct and exchange terms, in the DFT- 
LDA derived RPA it is given by the kernel 5h(r)/8p(v'). 

The system of equations [5] can be written in matrix 
notation as: 


= E^x[ k \ 

= (5) 



where we defined the RPA matrices A = {A^}, B = 
Ej-(E) = {Sij(E ±£j)} and the vectors X = {Xi}, 
Y = {Yi}. In the next subsection we will employ the 
ansatz of separability performing a transformation on 
this system of equations. 


B. Block-separable anzatz 

The standard separable ansatz states that for all 
ph,p'h! pairs, the two-body matrix elements can be fac¬ 
torized in terms of a one-body transition operator Q: 
A it j = A Q*Qj where Qi = (ip p \Q\<Ph)- Already in (tJ we 
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have relaxed this condition using three different coupling 
constants. We generalize this result to N blocks of ph 
excitations. In order to achieve such a description in an 
elegant manner, in the following we will adopt the matrix 
notation with compact indexes i = ph and j = p'hl . 

Let us consider the unique set of coupling constants 
A ij = Aij/(Q*Qj) = B it j/(QiQj). We construct the set 
of matrices and vectors: 


A — {A ij }, 

Q = { Qi j") 

z = ( Q* Q ) ( Y ) ’ 

x(E) = \Q\ 2 [£l 1 (E)-i+ 1 (E% 

where A is the real and symmetric matrix of coupling 
constants and Q the diagonal matrix of ph strengths. 
The Z vector hides more physical meaning being related 
to the transition amplitude for the k — th excited state 
(0|Q|fc) = |Z( fc )| = zj k ^ ■ The element of the diagonal 
matrix Xi(E) = 2|<5,;| 2 e i (£’ 2 — ef) can be identified as 
the linear response function of the system for a single ph 
state. 

Using the separable form of the interaction matrices 
A = Q*AQ and B = QAQ we can rewrite the system 
of equations [G] as: 

( Si 1 0 \ / Q*AQ QAQ \ /X\ = /X\ 
VO S+ 1 J V -Q*AQ* -QAQ* J V Y J V Y J ‘ 

Multiplying both sides with (Q*, Q) we obtain the sim¬ 
plified form: 


A = 


/ Aniii 

A 21121 


A 12112 

A 22122 


V Ajvii 


N 1 


• ' • AivllAf \ 
’ • • A2Afl2iV 

... Xin^nn / 


(9) 


where l nm is a matrix with all elements 1 and di¬ 
mension 7r n | x 1 7T m |, where 17r^ | is the cardinality of set 
7Tj. Thus, the system [7] is reduced to an equivalent 
A—dimensional system (with the same form), in which 
the quantities involved are now defined as sums over ph 
excitations in the same block Xn{E) = Yj^^XjiE), 

En = Yjeir n Eji IQImn = ^nmY2jen n \Q\j > n = ^> AT. 

The same system of equations can be obtained using 
the second quantization formalism and imposing the fol¬ 
lowing separability of the residual interaction: 


Sh 

Tp 

Zn 


1 

2 


'y ) XnmZnZrm 


— QjCj, 


( 10 ) 

( 11 ) 


where ct = a^dh- It can be shown that the vec¬ 
tor Z is related to the matrix elements of the op¬ 
erator z n by = (0|5 n |fe), and that the opera¬ 

tors z n verify the anti-commutation relation {ft, z m } = 
Snm{J2ie-K n IQ*I 2 ) = IQlnm Using the normalization of 
the initial RPA — |Y^| 2 ) = 1 and the de¬ 

pendency between Z and {X,Y} we can derive a nor¬ 
malization condition for the solutions of the [7] 


X(E) AZ = Z. (7) 

Note that this relation is equivalent with the original 
RPA equations [5] and since we defined the coupling con¬ 
stants for each ph^p'h 1 pair, no approximation has been 
done yet. Even though we halved the dimensionality of 
the problem, the new dispersion relation is more compli¬ 
cated, being contained in the condition of singular ma¬ 
trix: 

det (i - x(E) A) = 0. (8) 

Let us define the set of single particle excitations as 
7 r = {| ph)} and decompose it in N blocks n = (JiXi n i- 
We now employ the generalized concept of separability. 
We assume that pairs of particle-hole states interact with 
different coupling constant corresponding to the blocks 
they belong to Ay = A nm , V* G n n , j G 7r m , Vn, m = 1, N. 
With this ansatz one can rearrange the elements from [7] 
such that the A matrix becomes an array of constant 
blocks: 


Z^d E X~ 1 (E)Z = - 1 . ( 12 ) 

Recognizing that the elements of x are response func¬ 
tions for each block 7r n and having in mind the form of 
the residual interaction m we support the idea that the 
present formalism models a system of coupled subsystems 
of oscillators. 

Let us now stress out the main advantages of 
A-BSRPA: 

1. The dimension of the original problem 2 occ x unocc 
(which is infinite in principle) can be dramatically 
reduced to a system with dimension N. 

2. The number of solutions of the dispersion relation 
[5] is preserved since the latter can be written as 
an 2occ x unocc real polynomial in E. The rank- 
nullity prescribes one Z^ solution for each E^ 
excitation energy. 

3. The phenomenon of Landau damping is reproduced 
due to the presence of single particle effects in the 
block response function Xn{E). 
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4. The TKR (23j sum rule is recovered by the solutions 
of[7J 


m 1 {Q) = Y J \Qi\ 2 £i = Y.\ z{k) 

i k 


2 #. 


The goal of the whole RPA procedure is to reproduce 
the linear response function of the system. Since we pre¬ 
scribed Q as being the transition operator, the total re¬ 
sponse function (normalized to the sum rule) can be writ¬ 
ten as: 


| v -( fc )| 2 . . 

< 13 > 

where C(E,E^ k \T^) = 2E^[(E+iT^'>) 2 - E^ 2 ]- 1 
is the Lorentzian broadening of the delta function with 
some phenomenological width T^ -C E^. While this 
is the standard procedure in the representation of opti¬ 
cal spectra (for example) we should be aware that the 
block approximation involves basically an averaging over 
coupling constants, therefore, the broadening can be ex¬ 
plained as a feature of the dispersion of real \ p h,p'h'- 
We stress out that even if our results are very similar 
to the ones in [l] (Section II D), the physical picture 
is completely different: while they use a set of operators 
projected on the 7r space of ph excitations, we use a single 
operator projected on a set of n n spaces. Consequently, 
the interpretation of the coupling constants is different. 
The similarities arise from the mathematical apparatus 
which is basically the same. 


C. Particular cases 

Let us observe that the original separable RPA [3f is 
recovered if N = 1, i.e. one single block of relevant ex¬ 
citations is considered, which implies an unique coupling 
constant A between any pair of ph excitations. Then the 
dispersion relation [8] reduces to the well known equation: 


This means that for high energies, the interaction be¬ 
tween blocks of ph states can be neglected, in other 
words, the high energy excitations can be viewed 
schematically as collection of non-interacting subsystems 
(oscillators). These case can be justified keeping in 
mind that in clusters, the relevant single particle ener¬ 
gies (those with large strength \Q p h\) £ P h — leU while 
the collective states for the surface/volume plasmons are 
characterized by E a / V ~ SeV/^eV. 

In nuclear physics it is customary to consider the limit 
of a completely degenerate spectrum. This is motivated 
by the fact that the isotropic harmonic oscillator, which 
has equally distanced occupied levels of energy, repre¬ 
sents an empirically verified description of the mean field 
potential of the system in the ground state. Conse¬ 
quently, the relevant ph excitations are between equidis¬ 
tant shells. The same approximation is also valid in clus¬ 
ter physics. In this limit, £j = £,V* € 7r, the system of 
equations Q] becomes: 


|Q| 2 AZ = ^!_£! Z . (16) 

It follows that the solutions of the dispersion relation 
can be written as E^ = ±i/e 2 + 2 eqk, where {qk} are 
the eigenvalues of the | Q | 2 A matrix. Moreover, we obtain 
a simple normalization condition z[ : |Q| 2 Zfc = E^/eq^. 
Being obtained in the frame of degenerate limit, all N 
solutions describe collective excitations. 

At this end, we recognize that in our previous work 
we were not able to identify the connection between the 
analytical solutions derived in there and the fact that 
they were eigenvalues of a certain matrix (|Q| 2 A). 

Finally, we should link the generalized static polariz¬ 
ability for the Q operator with solutions provided by our 
formalism : 


N 


«(<9) = 

k 


IZ«I 2 

Ek 


(17) 


1 _ 2\Qi\ 2 £i 

A ~ ^ E 2 - £ 2 ' 

l 1 


(14) 


since it will serve as a supplementary test of the method. 


The N = 2 case has been treated in our previous work 
il]. At the other limit, when N = occ x unocc, i.e. each 
block contains a single ph excitation, the standard RPA 
is recovered. 

Also, the Tamm-Dancoff approximation can be ob¬ 
tained in our formalism by setting Y = 0 and defining the 
response function matrix as \ n = ^2 ie7rn \Qi\ 2 /(E — £j). 

We note that in the limit E £* the determinant 
from the dispersion relation [8] can be expanded in the 
first order and the equation becomes: 


III. RESULTS AND DISCUSSION 

Before presenting the numerical results, we will briefly 
describe the method used for the calculation of the 
ground state electronic spectrum. 

It is known Q that the response of clusters to exter¬ 
nal perturbations in the linear regimes is dictated by the 
valence electrons, while the core electrons and the nuclei 
remain frozen. The ionic background can be modelled as 
a positive charge distribution 


1 - Tr{x(E) A} = 1 - X] A nnXn(E) = 0. (15) 


n 


Pjeli r) = Po[l + exp(cr(|r| - R je i))} \ (18) 
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with 

Rjel = R() | L + ^2 J ■ (19) 

po is equal to the ratio 3n/4irr 3 , where r s is the Wigner- 
Seitzs radius. In practice, this value is adjusted such that 
Pj e i obeys the normalization condition f d 3 rpj e i(r) = n , 
where n is the number of atoms in cluster. The param¬ 
eters R 0 = r s n 1/,a and a™* = cq _m are a measure of the 
deformation of the cluster and cr replaces the diffusive 
character of the jellium which is known to give better 
results than the step profile (a —» oo). 

The ground state of the system was obtained by solving 
the DFT-LDA equations 

(^If v2 + v Ks(rj) <Mr) = (20) 

Vk-s(r) = V ion (r) + V H (r) + V xc {r). (21) 

Vion is the potential created by the ionic background 
(either from pseudopotentials, or the jellium model), 
the Hartree potential Vh is subject to Poisson equation 
V 2 Vjj = —47rp(r) and the exchange-correlation poten¬ 
tial V xc is approximated with the Gunnarson-Lundqviust 
functional 24]: 



r(Ro) 




Numerically, we solved the eigenvalue problems l20l for 
general clusters using a code that discretizes the equa¬ 
tions on a 3 D uniform grid with a step of ~ 0.11?o for a 
box with a length 5 times the radius Ro- The eigenprob- 
lem is solved using the imaginary time method coupled 
with the split-operator technique. For the kinetic op¬ 
erator, a Fourier representation has been employed, as 
well as for the computation of the Coulomb potential. 
For simplicity, we consider the jellium representation, ig¬ 
noring the ionic configuration. For more computational 
details see [25], [26|] . 

The residual two-body interaction prescribed bvEUlcan 
be written explicitly as: 


6h e 2 Sv xc (r) 

T, = yvTf + M7T 5(r - r > 


(23) 


Now, the numerical difficulty posed by RPA is obvi¬ 
ous, being reflected in the Ap^yh' terms which are com¬ 
puted as 6 D integrals for the Coulomb kernel. The LDA 
nature of the exchange-correlation potential makes the 
evaluation of matrix elements simpler as 3 D integrals. In 
principle, the long-range character of the Coulomb inter¬ 
action makes the assumption of separability impossible 
to prove. The level of applicability should be described 
only by numerical simulations of the true RPA matrix 
elements and comparison with their separable form. In 
practice the coupling constants A p h, P 'h' are computed for 
each ph — p'h' pair. 


FIG. 1: The electronic radial density, KS potential and 
first energetic eigenvalues are plotted for Na, 2 o and for 
singly charged Na^i as it follows: in the upper fig¬ 
ure the density solid line-A^a^j and dashed line Na 2 o', 
in the bottom figure the occupied states (red-iV 0 , 20 , 
orange-A r aJ 1 ) and the first 120 unoccupied states (blue- 
Nci 2 o, green- Na^) 


A. Sodium clusters 

The sodium clusters are known to be textbook cases 
in cluster physics due to their alkali character and for 
that we start with them as test cases. The value of 
r s « 3.93a.w. is used in the calculation of Rq. The 
spherical shaped clusters Na^ +1 ^ (with a magic num¬ 
ber .: n = 2,8, 20,40, 58,92,138,.. of atoms) are easier to 
tackle numerically, both from the perspective of DFT- 
LDA as well as of RPA calculations (for details, see (27|). 
Also, their optical spectra exhibits a strong surface plas- 
mon around 3eR 

The numerical challenge comes in deformed clusters 
where the natural degeneracies from the spherical case 
are removed and the spectrum of needed ph excitations 
considerably enlarged. This increases the dimension of 
the RPA system of equations and consequently the num¬ 
ber of matrix elements to be computed. In turn, as in 
nuclear physics, the deformation induces a splitting of 
the strong surface plasmon (landau fragmentation) 

As a demonstration of the numerical results regarding 
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the ground state we plot in Fig. |T]the radial profile of 
the electronic density in Na 2 o and Na^ as well as the 
first few energetic levels. 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 


E(eV) 

FIG. 2: Single particle ph excitation strengths in Na §g 
for the dipole operator (blue) and volume operator 
(red)) 

The formalism described in the previous section is gen¬ 
eral and valid for any transition operator Q as long as 
the block-separable anzatz holds. In practice we should 
use the standard dipole operator to describe the optical 
spectrum, or, a more complex local operator of the fol¬ 
lowing form: Q{ r) = TO r z l) m (n) j28j. This operator 
accounts for dipolar surface oscillations and so, for the 
surface plasmon. To account for the volume plasmons 
we should employ operators which have non-zero lapla- 
cian character, generically Q(r ) = r p Yi m (fl),p ^ l. In 
the end we use the operators Q s = z for dipole normal 
modes and Q v = r 2 for volume modes.!!! 

Information about relevant excitations in the system 
can be obtained by computing the spectrum S(E) for the 
independent particle response. In the limit of A —>• 0, 
the spectrum [l3] reduces to the trace of In Fig. [2] 
we have plotted this quantity for Nci^i both for Q s and 
Q v . The same behaviour can be seen in all sodium clus¬ 
ters. The degenerate limit is justified because of the nar¬ 
row spectral region in which the matrix transition oper¬ 
ator exhibits large values. In particular we can com¬ 
pute the degenerate ph energies as a weighted mean: 

S s/v = J2ph £ ph\Q s p / h v \ 2 / T, ph \Qph?- 

In order to prove the validity of the separable ansatz, 
we perform an analysis of the A p h, P 'h’ set in respect with 
the characteristics of ph states. Since the obtained val¬ 
ues span many orders of magnitude we have plotted in 
Fig. [3] the logarithmic dependence between A p h,p'h' and 
\Q p h\,\Q P 'h'\ f° r Nag. Surprisingly enough, we have 
found that the set of data points lies roughly in a plane. 
Moreover this type of behavior has been found in all the 
other simulated clusters. This has driven us to propose, 
without any mathematical justification, the following em¬ 
pirical relation A p h, P 'h' = ciQphQp'h') 1 , where c > 0 and 

7 < 0 . 



^9 (I Qph I) q h'\) 


FIG. 3: Scattered plot in Nag for lg(\X p h,p'h'\) vs. 
lg(\Qph\),lg(.\Qp’h'\) fitted with a plane 


With these, it becomes obvious that the validity of the 
block ansatz relies on the possibility of defining the blocks 
from the distribution of Q p h elements. This is visually 
shown in Ufor Na£ in a logarithmic histogram of \Q p h\ 
where the a choice of N = 4 blocks is represented in a 
colour plot. 

As a conclusion of these features provided by the nu¬ 
merical computation of RPA A p h, P 'h' coupling strengths 
we propose the following scheme: 

• Solve the KS problem for the system and obtain 
the occupied and unoccupied spectrum. 

• Compute the {s p h,Q P h} strengths for the single 
particle ph excitations. 

• Choose the ph blocks Tr n separating regions in the 
distribution of Q p h strengths. 

• Using three free parameters c, 7 compute the block 
coupling constants X nm = c((Q) n (Q) m ) 1 , where 
{Q) n = Sie-n-n Qi/Knl is the mean in the block. 

• Solve the dispersion relation [5] and the system of 
eqns. [7] for block configuration 

• Compute the polarizability accordingly with [T7] 

• Tune the free parameters such that the resulting 
polarizability ct(Q) approaches the experimental 
value. 

• Compute the response function from 1131 

In order to appreciate the correct magnitude of the 
coupling constants, we first apply the A’-BSRPA in the 
degenerate limit of the case N = 1. The degenerated 
energy is computed as above e = mi(Q)/mo(Q). We take 
Q = ez the dipole operator and using the solutions for 
IT7)1 eigenvalue problem we compute the A that reproduce 
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FIG. 4: Distribution of /glQpfcJ in Na g" for the dipole 
operator. The possible configuration of 4 blocks 77 , * = 
1,4 is represented with different colours. 
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FIG. 5: Values of A in the case of one-block and degen¬ 
erate limit obtained from the DFT-LDA prescription 
for the degenerate energy and the experimental values 
for the peak of the surface plasmon [3(| 



both the static polarizabilities (taken from j 2 E|) and the 
centroid of the surface plasmons. The clusters Na n with 
n = 8 , 20,40, 58, 92,138,196 have been chosen due to the 
fact that in deformed ones, the degeneracy is not a valid 
approximation anymore. The computation reads A = 
(E^ — e 2 )/(2emo) and the results can be seen in Fig. [5] 
From a numerical fit, we have found the following values 
for the parameters of the empirical formula 7 ss —0.75 
and c « 0.025eV/(e 2 Oo) (where e is the electron’s charge, 
ao the atomic unit of length and Q p h is expressed in ea 0 ). 

Going a step further, we disregard the degenerate limit 
and solve the 1-BSRPA problem for the same neutral 
spherical clusters obtaining the optical spectra plot in 
Fig. [ 6 ] in respect with the true RPA results. A generic 
width Tfe = 0.02eV has been used to broaden the reso¬ 
nances accordingly with eq. [13] as the lorentzian param¬ 
eter. Also, we have used for A the values obtained in Fig. 
[5] From the numerical results, a gross agreement can 
be found between full RPA and V-BSRPA even for this 





FIG. 6 : Optical spectra obtained with full RPA (red, 
dashed) and 1-BSRPA (blue) in the one-block case for 
spherical clusters Na n where n = 8 , 20,40, 58,92,138 


minimal case N = 1. The differences can be explained 
recognizing that N = 1 case is a crude approximation. 

In order to have a final test for 1V-BSPRA, we in¬ 
vestigate the optical response for singly charged sodium 
clusters Va+, n = 9,21,41,59 which are also spherical 
symmetric. Even though, in principle, the most accu¬ 
rate results should be recovered for N -4 occ x unocc, 
we have found that the formalism with N = 4 provides 
the best computational cost-accuracy ratio. For that rea¬ 
son we compute in Fig. |7]5(E) for the above mentioned 
clusters employing N = 4 and using the discussed em¬ 
pirical parametrization for X nm . The parameter c has 
been left free and tuned in order to obtain the RPA pre¬ 
scribed static polarizability as a first test. The results are 
compared with full RPA calculations which are in good 
agreement with the experimental data j3lj. We conclude 
at this end that V-BSRPA gives a pertinent description 
of the dipole excitations in sodium clusters. 


B. C60 Fullerene 


Since the discovery of Cq 0 fullerene [32] a lot of ex¬ 
perimental and theoretical attention has been drawn to¬ 
wards its spectacular properties. In particular, its optical 
spectrum is known to be dominated by two surface plas¬ 
mons around 20eV and 40eV (33}. In electron energy 
loss experiments three volumes plasmons emerge in the 
20 — 30eV range f3l| . The theoretical work done in this 
direction has been able to reproduce the experimental 
values from many perspectives as RPA [35], ab-initio [], 
TDDFT [33[, Time-Dependent Thomas Fermi |36j], etc. 



































FIG. 7: Optical spectra obtained with full RPA (red, 
dashed) and with 4-BSRPA in Nat, where n = 
9,21,41,59 


Since it is genuinely a more complicated system than 
the simple alkali Na clusters we considered it as a good 
further test case for the present formalism. 


In order to solve the KS equations for the ground 
state we have modelled the ionic background accord¬ 
ingly with [H, [37]] defining a jellium shell pj e i( r) oc 
@(r — ri) 0 (r 2 — r) centred on the position of the car¬ 


bon nuclei, where rq = 2.8 A, r 2 = 4.2 A and normalized 
to n = 240 which is the number of considered valence 
electrons. Following j38| we have added a supplementary 
pseudopotential V 0 to the total KS potential in order to 
account for the ionic structure. Due to the spherical sym¬ 
metry of the cluster of Cqq we were able to solve numer¬ 
ically only the radial KS problem and so, a more refined 
grid. For the calculation of matrix elements we have used 
the multipolar expansion of the Coulomb kernel (for more 
details about calculations of RPA in spherical systems see 

13 ). 


A first set of results is plotted in Fig. [5] where the 
electronic density, the effective KS potential and the first 
few energetic levels are represented. 


Using the prescription for A obtained for sodium clus¬ 
ters we use the 1 — BSRPA for Cgo to construct the op¬ 
tical spectrum in the case of dipole operator. The result 
can be seen in Fig. ?? where the main surface plasmon 
from around 20eU is well reproduced. A spurious split¬ 
ting is present as well as a second, smaller peak around 
32eU which can be interpreted as the surface plasmon 
from around 38eU. 


Finally, we use the 5-BSRPA with the Q s (j) = ez 
for the optical spectrum. The same parametrization has 
been used for A ram as in the case of sodium clusters, and 
the free parameter c has been tuned to reproduce the 
experimental static polarizability f3f|. In Fig. [TUlare re¬ 
produced the numerical results vs the experimental spec¬ 
trum [4p| . 



r(Ro) 



FIG. 8: In figure a) is plotted the density of electrons 
(solid line) vs. the density of jellium (dashed line); Fig¬ 
ure b) presents the KS potential (solid line), the occu¬ 
pied (red lines) states and first few unoccupied (blue 
lines) energetic levels 



E(eV) 


FIG. 9: Dipole spectrum for C 60 obtained for IV = 1 
block and A ss 0.019 


Conclusions 


Starting from DFT derived RPA equations we employ 
the ansatz of separable residual interaction defining N 
blocks of ph excitation with specific coupling constants 
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between them. We show that, in the frame of this approx¬ 
imation, the RPA eigenvalue problem can be converted 
in a homogeneous system of N equations. Associated dis¬ 
persion equation and normalization condition is derived 


for the unknowns of the reformulated problem. 

The present WBSRPA formalism has the main advan¬ 
tage over the standard RPA that reduces dramatically 
the dimensionality of the eigenvalue problem conserving 
the number of the solutions and consequently the Landau 
damping feature of the normal modes. Also, the nature 
of the collective excitations is naturally reflected in the 
degenerate limit of the system. On the other side, the for¬ 
malism does not prescribe a recipe for the computation 
of coupling constants, which have to be tuned in order 
to reproduce the experimental data. Still, an empirical 
formula is proposed which holds in the case of sodium 
clusters and Ceo fullerene, but we do not sustain the idea 
that the parametrization is universal. 

We conclude that the present generalization of the sep¬ 
arable RPA is a real alternative to the high costs of 
the full RPA calculations, but the problem of coupling 
constants remains open. We support the idea that the 
present formalism could find applications in many other 
systems from nuclear to solid state physics, but also in 
modelling other phenomena as the pairing is. 
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